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The Gaussian Approximation, proposed originally by Ottinger [J. Chem. Phys., 90 (1) : 463- 
473, 1989] to account for the influence of fluctuations in hydrodynamic interactions in Rouse 
chains, is adapted here to derive a new mean-field approximation for the FENE spring force. This 
FENE-PG force law approximately accounts for spring-force fluctuations, which are neglected 
in the widely used FENE-P approximation. The Gaussian Approximation for hydrodynamic 
interactions is combined with the FENE-P and FENE-PG spring force approximations to obtain 
approximate models for finitely-extensible bead-spring chains with hydrodynamic interactions. 
The closed set of ODE's governing the evolution of the second-moments of the configurational 
probability distribution in the approximate models are used to generate predictions of rheological 
properties in steady and unsteady shear and uniaxial extensional flows, which are found to be 
in good agreement with the exact results obtained with Brownian dynamics simulations. In 
particular, predictions of coil-stretch hysteresis are in quantitative agreement with simulations' 
results. Additional simplifying diagonalization-of-normal-modos assumptions are found to lead to 
considerable savings in computation time, without significant loss in accuracy. 
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I. INTRODUCTION 

The importance of hydrodynamic interactions (HI) in 
determining the behaviour of dilute polymer solutions 
has been widely recognized since the introduction of the 
Zimm theoryi. The limitation, however, of the accuracy 
of Zimm theory to linear viscoclastic property predictions 
has spurred the development of more sophisticated treat- 
ments of HI. For instance, while the consistent averaging 
approximation^ replaces the equilibrium averaged Oseen 
tensor in Zimm theory with a non-equilibrium average, 
the Gaussian approximation^ further accounts for fluc- 
tuations in HI. The Gaussian approximation, which is 
the most sophisticated approximation to date, has been 
shown to be highly accurate in shear flows by comparing 
its predictions of viscometric functions with exact Brow- 
nian dynamics simulations (BDS)^. 

The use of bead-spring chain models with linear 
Hookean springs in all these models has prevented the 
evaluation of their accuracy in extensional flows, in which 
the polymer coil undergoes a coil-stretch transition at a 
critical value Wic of the Weissenbcrg number Wi = Xe, 
where A is a characteristic relaxation time, and e is the 
extension rate. As is well known, the extensional vis- 
cosity becomes unbounded at Wic because of the infinite 
extensibility of Hookean springs. The prediction of a 
bounded extensional viscosity was achieved in early the- 
ories by replacing Hookean springs with nonlinear finitely 
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extensible springai^. However, these theories did not in- 
corporate HI. The important role played by HI in exten- 
sional flows has been established beyond doubt in a series 
of recent publications, which have shown that the inclu- 
sion of HI leads to quantitatively accurate predictions 
of both macroscopic and mesoscopic properties^'S'^iiSiii. 
All these theoretical predictions have been obtained by 
carrying out exact BDS incorporating fluctuating HI and 
nonlinear finitely extensible springs. 

An additional aspect that has been shown in these 
works to be crucial to obtaining quantitative predictions 
is the need to use bead-spring chain models with a sufR- 
ciently large number of springs, Ng (or equivalently, de- 
grees of freedom). This requirement is eventuated by the 
need to accurately reflect both the changing nature of 
the drag experienced by the polymer as it unravels in an 
extensional flow, and the self-similar character of a high 
molecular weight polymer (which is responsible for the 
observation of universal behaviour in dilute polymer so- 
lutions). The requirement of large Ng, coupled with the 
need to include HI (which leads to the CPU time scaling 
as Ng-^), makes the use of BDS for routine calculations 
highly impractical. While the development of schemes 
for the accelerated calculation of HI has alleviated this 
problem somewhati^, there is a pressing need for the de- 
velopment of approximations capable of describing exten- 
sional flows, that are both accurate and computationally 
inexpensive. The aim of this work is to introduce such 
an approximation by extending previous approximations 
for HI, which were based on bead-spring chain models 
with Hookean springs, to models with nonlinear flnitely 
extensible springs. 

Wedgewood and Ottinger— have previously introduced 
a similar approximation by combining the consistent av- 
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eraging scheme for HI with a FENE-P approximation for 
the springs. While they carried out a detailed investiga- 
tion of the predictions of the model in shear flow, they 
did not evaluate its accuracy by comparison with exact 
BDS, nor did they obtain any predictions in extensional 
flows. In this work, we introduce an approximation that 
extends the previously established Gaussian approxima- 
tion by accounting for fluctuations in both HI and the 
spring forces. The accuracy of this model is then verified 
by comparison with BDS in a number of both transient 
and steady, shear and extensional flows. 

The reduction in computational time achieved by the 
introduction of such approximations, despite being signif- 
icant, is not sufficient to explore the behaviour of truly 
long chains necessary to adequately model high molec- 
ular weight polymers. Several workers have shown that 
mapping the bead-connector vectors to 'normal' coordi- 
nates, followed by an assumption that certain matrices in 
the transformed space are diagonal, leads to further sig- 
nificant reductions in CPU time, without any sacrifice in 
accuracjii^ii^. Here, we carry out a similar diagonaliza- 
tion procedure, and examine both the reduction in CPU 
time, and its effect on the accuracy of the model. 

The paper is organized as follows. In the follow- 
ing section, we recall the basic equations for non- 
free-draining bead-spring models of dilute solutions of 
finitely-extensible polymer molecules. Subsequently, we 
develop in section IIIII the equations for the Gaussian 
approximation, and its diagonalized version. Section If VI 
introduces the flows and material functions of interest 
in this paper. The symmetries in the flow patterns con- 
sidered lead to considerable simpliflcations in the imple- 
mentation of the numerical methods for computing the 
material functions, which are briefly discussed in section 
Ivl Predictions of the approximations for the material 
functions in shear and extensional flows are presented in 
scction lvTl along with the exact results obtained with BD 
simulations. Section IVIII summarizes the central conclu- 
sions of this study. 



II. BASIC EQUATIONS 

Under homogeneous conditions, the configurational 
state of a bead-spring chain, with Ng springs, is com- 
pletely specified by the set of connector vectors {Qi \ i = 
1, . . . , Ns]^. The Fokker-Planck equation governing the 
evolution of the configurational probability distribution 

i'iQi, - ■ ■ ,QN,,t) is. 
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where, is Boltzmann's constant, k is the transpose of 
the position- independent traceless velocity gradient \7v, 



T is the absolute temperature of the solution, and C is 
the Stokesian bead- friction coefficient, which is related to 
the bead radius a through ( = Swarfs, with ijs being the 
viscosity of the Newtonian solvent. The dimensionless 
diffusion tensors Aij are discussed in detail shortly. The 
spring connector force F^' is given by. 
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with Sj representing the potential energy of spring j. 
Here, the following FENE expression for the spring force 
is usedii: 
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Qj=HaQj)Qj, (3) 



where, H is the spring constant, Qj = \Qj\, Qo is 
the 'fully-stretched' length a single spring, and ^ = 

[l - Q] /Qo] " is the nonlinearity in the FENE spring 
force law. 

By considering a bead-spring chain of A's springs to 
represent a polymer molecule of TVk Kuhn-segments, each 
spring can be imagined to represent a "sub-molecule" of 
-^K. s = N-n/Ng Kuhn segments. The parameters in the 
FENE spring- force law, H and Qo, can then be related to 
the two parameters characterizing the equilibrium struc- 
ture of a sub-molecule, namely, &k (the length of a Kuhn 
step) and A'k, s- Since the contour length of the sub- 
molecule is BkNk, s, it follows that Qo = 6k-/Vk. s. The 
equilibrium mean square end-to-end distance of a sub- 
molecule, in the absence of excluded volume forces, is 
i?g = ^i-^K, s- Using the equilibrium Boltzmann distri- 
bution for a single spring governed by the FENE spring 
force law, tpl^ ~ exp(— S'^/ZcbT), it can be shown tha1ii&, 
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Substituting for and Qq in terms of and A^k. s, leads 
to the following expression for the spring contant H, 
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The equations above indicate that taking the limit 
Qo ^ oo, while keeping i?| constant, leads to ^ 1. 
In this fimit, Rl = 2,k^T/H, and the FENE spring- 
force expression reduces to the Hookean spring force 
law, F^^'^ = HQ. In this work, chains with Hookean 
springs are referred to as "Rouse" chains, while finitely- 
extensible bead-spring chains are denoted as "FEBS" 
chains. ^ 

The diffusion tensors Aij in Eq. are given by. 



(6) 



where, A 
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<^|i-j|,i is an element of the A^s x Rouse matrix, and 



3 



^u^i, i^, ^ = 1, . . . , iV, arc HI tensors, which are typically 
related to the inter-bead displacement r^^ between beads 
fi and ly, by expressions of the form, 
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Note that the inter-bead displacement vector r^j, can be 
related to the spring connector vectors through. 



sgn(i^ - pi) 
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Models which neglect the influence of HI — such as 
the original bead-spring model of Rousei^ — can be inter- 
preted as ones in which the tensor C is set to zero. Early 
models incorporating HiiiiS used the Oseen-Burgers form 
of the HI tensor, in which 
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Although the Oseen-Burgers description of HI is inac- 
curate when the pair-wise separation between beads is 
comparable to the size of the beads, the hydrodynamic 
behaviour of long polymer chains is dominated by the 
interactions between non-adjacent beads, and the short- 
range inaccuracy of the Oseen-Burgers description has 
been shown to be relatively unimportant in the predic- 
tion of properties that depend on the polymer chain as a 
wholaSiSa. 

Further, its simple form is particularly useful for the 
development of closure approximations. In numerical 
simulations, however, the Oseen-Burgers expression be- 
comes problematic when r^^ < a, since the 3N x 3iV 
diffusion block-matrix comprising the Aij tensors as its 
constituent blocks, can become non-positive dcfinitoSi 
when beads overlap. In its place, in all the Brownian 
dynamics simulations performed here, the Rotne-Pragcr- 
Yamakawa (RPY) modificationSiiS^, 
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rn^ > 2a, 
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which ensures that the diffusion block-matrix always re- 
mains positive-definite, is used. 

In order to express the Fokker-Planck equation in di- 
mensionless form, the following length scale and associ- 
ated time scale are introducedSiiS, 



Ik 
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In the Hookean limit, 
conventional form, 



4fcBr' 

both £s and Ag reduce to their 
= ia = y/k^T/H and As = Ah = 



C/4ff = Cel/Ak^T, respectively. Defining t* = t/As, 
Q* = Qi/ £s, K* = AgK, the Fokker-Planck equation can 
be recast in dimensionless terms as 



dt* 
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where, 



H* = 
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[12) 



(13) 



is the dimensionless spring constant. When €h is used 
as the length scale, the parameter Qg is represented in 
dimensionless terms by the finite extensibility param- 
eter, h = Q\l IS^. The counterpart of this parame- 
ter in the non-dimensionalization scheme adopted here 
is, h* = Qo/is- Since Qo = 6k^k, s, it follows that, 
b* = 3iVK, s- Because of the direct relation of b* to the de- 
gree of coarse-graining A^k. s, the latter is often reported 
instead of b* when presenting results. 

On introducing the hydrodynamic interaction param- 
eter, h*, defined by. 



(14) 



the dimensionless tensor C^i^fj, occurring in the defini- 
tions of the diffusion tensors [Eq. ©] can be written as 



30F/1* 



4r 



(15) 



where r*^, = (1/ £s)rf,u- 

All the rheological properties of interest in the present 
work-which are discussed in greater detail in Section llVl - 
can be obtained from Kramer's expression for the poly- 
mer contribution to the stress tensor T^p^, which in non- 
dimensional terms is, 



knT 



rp^N,S-H*J2{i,Q*Q* 



(16) 



Here, rip is the number density of polymers. 6 is the unit 
tenor, and angular brackets denote an average performed 
with the configurational distribution function ijj. 

Besides rheological properties, the mean square end-to- 
end distance R* ^ is used here to provide information on 
the size and shape of polymer molecules as they stretch 
and orient in response to the imposition of flow. 
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Further, the prediction of the various models of the bire- 
fringence in extensional flows is evaluated using the fol- 
lowing expression^Si^: 



a multivariate Gaussian, 



An* 



1 



{Q*l 



(18) 



where, n' is the intrinsic birefringence tensor, and A is a 
constant that, besides depending on the chemical struc- 
ture of the polymer and the solvent, is proportional to 
the polymer concentration Up. 

Eauations llt)lll7l andll8labovc show that all the macro- 
scopic properties of interest here can be described in 
terms of expectations of configuration dependent func- 
tions. It is possible to obtain formally the equation for 
the time evolution of any configuration dependent func- 
tion by multiplying the Fokker-Planck equation with the 
function and then integrating over all possible configura- 
tions. However, as a direct consequence of the nonlinear- 
ities in the model arising from FE and HI, the resulting 
equations are not closed with respect to the expectations 
whose evolution equations are desired. For instance, the 
following set of evolution equations for the second mo- 
ments of the probability distribution can be derived in 
this manner: 



: q:q^ )(i) = ^^^^ Q*Q*3 ) - ( ) • X* 



H* 



(19) 



where the subscript "(1)" denotes the contravariant con- 
vected time derivative of a tensor—. As a result of the 
nonlinearities in the functions ^ and O introduced by 
the spring force law and HI respectively, the equation 
above for {Q*Q* ) is not closed with respect to the sec- 
ond moments. It is this 'closure problem' that has moti- 
vated both the introduction of various closure schemes, 
and the development of exact Brownian dynamics simu- 
lations for obtaining predictions of properties away from 
equilibrium. As mentioned previously, the Gaussian ap- 
proximation for Rouse chains is the most accurate closure 
approximation to date. In the section below, we develop 
the Gaussian approximation for FEES chains with HI. 



III. GAUSSIAN APPROXIMATION 

For Rouse chains, Ottinger— and Wedgewood^ pointed 
out that the complicated averages on the right-hand side 
of equation 1191 above (with = 1) can be expressed as 
functions of the second- moments {Q*Q*) by assuming 
a priori that the configurational distribution function is 



ip = J\f exp 
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Q* 
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(20) 



where A/" is a normalization factor. The tensors p*^ are 



related to the second moments a* 



[q:q; 



of the 



Gaussian distribution through the following set of 
linear algebraic equations 



E 



= SijS 



(21) 



As is well known, this implies that the Gaussian distri- 
bution -0 is completely specified once its second moments 
cr*^- arc known. 

The crucial property of Gaussian distributions that en- 
ables the simplification of complex moments is described 
by Wick's decomposition rule. In the present instance, 
because the right hand side of Eq. H19|l has complex av- 
erages of the form {xxf{x)), where x is the random vari- 
able of interest, and / is a nonlinear function of x, Wick's 
decomposition rule leads to 



df 

{xxf) = {xx){f) + {xx) ■ {-g^x) 



(22) 



By applying Eq. to Rouse chains with HI, lOttingeJ 
showed that Eq. ifT^ (with ^„ = 1) reduces to a closed 
set of equations for the second-moments {Q*Q*)- In 
principle, the same approach can be used to obtain clo- 
sure of Eq. H19|) with both the nonlinear phenomena 
of finite extensibility and hydrodynamic interactions in- 
cluded. However, directly using the Gaussian assumption 
and Wick's decomposition to simplify the expectations on 
the right-hand side of Eq. (|19|) leads to the incorrect pre- 
diction of a dependence of equilibrium static properties 
on h* . Similar behaviour was earlier observed by Prab- 
hakar and Prakash— when the Gaussian approximation 
was used to achieve closure in a Hookean dumbbell model 
simultaneously incorporating excluded volume and hy- 
drodynamic interactions. The reason for this incorrect 
prediction is that, on direct application of the Gaussian 
approximation, the contributions from FE (and/or ex- 
cluded volume) and those from HI cannot be separated 
into independent factors at equilibrium—, thus leading 
to an unphysical dependence of the predicted equilib- 
rium second moments on the HI parameter, h* . To avoid 
this problem, closure of Eq. H19|) is achieved here in two 
steps. In the first step, a mean-field quadratic spring 
potential 5^, which approximately accounts for fluctua- 
tions in the FENE spring force, is derived by using the 
Gaussian approximation for free-draining FEBS chains. 
A bead-spring chain model is then constructed, which 
incorporates hydrodynamic interactions, and which has 
springs described by S'„j. In the second step, the Gaus- 
sian approximation is again invoked to obtain a closed 
set of equations for the second moments of this model. 
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In this manner, the incorrect prediction of a dependence 
of equilibrium static properties on h* is avoided. 

Before discussing these two steps in detail in sec- 
tions and below, it is appropriate to note that 
if / had been replaced with its average (/) in Eq. 1221 
thus neglecting the fluctuations in /, one obtains the first 
term on the right-hand side alone. Thus, for a Gaussian 
distribution, it is seen that fluctuations in / are exactly 
accounted for by the second term on the right-hand side 
of the equation above. This convenient separation of the 
complicated expectation on the left hand side into contri- 
butions from the mean and fluctuations will be used here 
to highlight the important role played by fluctuations in 
the spring force and in hydrodynamic interactions. 



A. Free-draining FEBS chains 

In this section, the Gaussian approximation is used 
to achieve closure of the second-moment equations for 
chains with finitely-extensible springs, but without HI. 
Prakash and co-workers have previously used the Gaus- 
sian approximation for bead-spring chain models without 
HI, but which are, however, still non-linear because of 
the presence of conservative non-hydrodynamic excluded 
volume forces between beads. The approximation was 
found by them to lead to qualitatively accurate prcdic- 

tionol5,28i29,30. 

For free-draining chains, Aij = AijS, and hence from 
Eq. H19|l. one obtains for FEBS chains 



in place of the exact spring potential S*^. In other words, 
applying the Gaussian approximation to simplify the av- 
erage due to FENE springs in free-draining chains is 
equivalent to using the quadratic mean-field spring po- 
tential S'„j, since this approximation leads to 
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(23) 



The equation for the polymer stress is the same as in 
Eq. (|16|l . In this case, closure of the evolution equation 
above is not possible because of the FENE nonlinear- 
ity, ^rrn On the right-hand side. However, as in the case 
of the Gaussian approximation for HI, assuming that 
the configurational distribution is a Gaussian, and us- 
ing Wicks decomposition rule, the problematic terms on 
the right-hand side of the equation above can be for- 
mally simplified. Since dS^/dQ^ = H*(^rnQ*m (where, 
S*^ = S,nlk^T), using Eq. ^ leads to 



q: 



(24) 



The expression on the right-hand side could have also 
been obtained by directly using a quadratic mean-field 
spring potential. 



S = 
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Prakashi^ have previously shown that Fixman's treat- 
ment of the excluded volume interaction potential^ cor- 
responds to the use of a similar pair-wise mean-field 
quadratic excluded volume potential. 

The expression for in Eq. [251 can be simplified by 
using dS;jdQ*„, = H*UQ*m. 



/ di„ 

\dQ; 



-q: 



Q*mQ 71 



(27) 



The first term within the brackets on the right-hand side 
of the equation above containing {£,m), the average of the 
spring force nonlinearity, can be recognized as that aris- 
ing from a consistent-averaging treatment of that nonlin- 
earity, while the remaining term accounts for the influ- 
ence of fluctuations in the spring force. 

Recursive application of Wick's decomposition rule, us- 
ing the result 9Cm/<9Q*„ = 2Q*„ d£_,n/d{Qm'^), leads to a 
series expansion, 



H* 



(6> 



>6 + ^2^ 



d{Q* 



0-: 



• QmQn 



(28) 



where 0*^^ 



denotes the matrix product of s copies of 
Omrn- The expansion above is valid provided the Gaus- 
sian averages (^m) and {d^^ml d{Q*my) exist for all 
s = I, . . . , oo. For FENE springs, however, 

I 



(1 



)* 2 



and. 



d{Q*ly b*^ {i-Q*„yb*) 



\s+l 



(29) 



(30) 



The Gaussian averages of these functions do not ex- 
ist since all of them have non-integrable singularities at 
Q*n — whereas the Gaussian distribution has an in- 
finite range. To make analytical progress, it is necessary 
to introduce further approximations. Here, analogous to 
the Pcterlin closure approximation^. 



(^m) is replaced with 



1 



{{l-Q-LVb*)) 



(31) 



and. 



" Sri 



diQ* 



is replaced with 



1 



(32) 
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Although each of these terms may remain bounded as 
long as (Q^^) = tra*„„ < 6*, there is a need to care- 
fully analyze the convergence properties of the infinite 
series in Eq. H28|) after introducing these approximations. 
Here, only the first term in the infinite series in Eq. H28() 
is retained to represent spring force fluctuations. As a re- 
sult, the following mean- field quadratic spring potential 



The use of this potential leads to a reduction of the com- 
plex moments in Eq. 1231 

H*{Q*QUm)^H*al^-L^, (35) 

and a mean-field spring force given by, 

- H*L„, ■ Q;, . (36) 

As mentioned above, neglecting the spring force fluctu- 
ation term in i,„ is equivalent to a consistent-averaging 
treatment of the FENE nonlinearity, which leads to the 
well known FENE-P force law proposed originally by 
Bird et alM, 

= ^ \-frl* //.* ^"' = H%,Q:n (37) 

where, Cm = [1 — tr Incorporation of fluctu- 

ations through the Gaussian approximation leads to the 
spring force expression in Eq. H3t)|l . denoted here as the 
"FENE-PG" spring force law. 

The evolution equation for the second moments in the 
Gaussian approximation for free-draining FEBS chains is 
consequently 

^ij, (1) ~ ^ ^ i'^im ■ LmAmj + AimLjn. ' tT*y) 

m—1 

(38) 

and the polymer stress predicted by the approximation 
is given by 

rl^N,b~H*Y,al-L,. (39) 

1=1 

It is worth noting that for finite cr^^, as 6* ^ oo, the 
fluctuation contribution in Eq. (|34|1 becomes vanishingly 



is used, 

^m = ~:^I^m' QmQm i (33) 

where. 



(34) 

I 

small, Lm — 6 and the spring becomes more Hookean- 
like. On the other hand, as the bead-spring chain ap- 
proaches its full extension, and 

fb* 0\ 
<m ^ 0, (40) 
\0 oj 

the scalar prefactor in the fluctuation term in Eq. I|34() 
approaches (1/6*), which however becomes negligible 
compared to the FENE-P contribution, which diverges. 
Thus, as chains become highly stretched, predictions ob- 
tained with the FENE-P or FENE-PG approximations 
are expected to be nearly identical. 

Although the FENE-PG approximation is new and 
has been developed in this study, the FENE-P approx- 
imation has been studied in great detail by several 
workera^i^'^i'^"^i^^i^^i^'''i'^^i^^. As pointed out in many of 
these studies, perhaps the single most important dif- 
ference between the exact FENE and the mean-field 
FENE-P and FENE-PG approximations lies in the na- 
ture of the configurational probability distribution func- 
tion. The singularity in the non-averaged FENE spring 
force strictly restricts the lengths of the connector vectors 
to lie in the domain [0, Qo)- However, with the Pcterlin 
closure, the solution of the modified Fokker-Planck equa- 
tion with FENE-P/PG springs is a Gaussian. Therefore, 
there is a non-zero probability that a spring has a length 
greater than Qq. This has been clearly illustrated by 
Keunings^ who used Brownian dynamics simulations of 
dumbbells with FENE-P springs to show that in strong 
cxtcnsional flows, a significant proportion of the dumb- 
bells had length greater than the "maximum allowed" 
dimcnsionless extension Qq. However, in all the litera- 
ture on the FENE-P model, the mean-squared lengths of 
the springs are observed to always be lesser than Qq. In 
other words, the mean-field FENE-P model appears to 
satisfy the maximum extension constraint in an average 
sense. 

It is observed in the studies cited above that predic- 
tions with the FENE-P approximation compare well with 
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1 - tr o-*,„„/6* 1 - 2 tr a;„Jfe* + [( tr a*„^ J^ + 2a*,„„ : a*„ J/b* ^ ' 
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the FENE model for steady-state properties, in both 
shear and uniaxial extensional flows. An explanation for 
the good agreement at steady-state in a strong flows is 
based on the observation that under such conditions most 
of the springs in the chain are highly stretched. The 
exact steady-state distribution function for the length 
of the spring in a dumbbell model is sharply peaked 
around ^ (Q^), and can be well approximated by a (5- 
function^iiS. However, during start-up of shear and ex- 
tensional flows, predictions with the FENE-P model for 
the transient rheological properties deviate strongly from 
the exact results of BD simulations prior to the attain- 
ment of the steady-state2£*S. 



of algebraic equations: 



m—1 



1 



(41) 



At equilibrium, since cr*^ — a*^ ^^3 and L„i, cq = 
Lm. cq3, the equation above can be rearranged into the 
following set of equations. 



(42) 



There have been several previous attempts to improve 
the FENE-P approximation. In conditions where springs 
can be highly stretched, the true spring length distri- 
bution is better approximated as the sum of a uniform 
distribution for smaller values of Q and a (5-function lo- 
cated at a large Q closer to (but less than) Qo- This 
FENE-L approximation, and its three-dimensional coun- 
terpart, are shown to be more accurate than the FENE- 
P model, particularly in predicting the phenomenon of 
stress-conformational hysteresia^Miii^. Grban et alm^ 
show that an entropic argument can be used to im- 
prove upon the FENE-P approximation by first expand- 
ing the exact FENE potential in a Taylor's series ex- 
pansion about (Q^) and then systematically regrouping 
terms in the series before truncating it. This approach 
gives rise to a hierarchy of spring potentials of increas- 
ing complexity, the simplest of which is shown to be the 
quadratic potential corresponding to the FENE-P force 
law. By using the next higher order potential (which 
is denoted as the "FENE-P-l-1" potential by Gorban et 
al. ) in the hierarchy, these authors demonstrate with 
a one-dimensional dumbbell model that accuracy is sig- 
nificantly improved in an unsteady elongational flow. In 
Section IVll below, the accuracy of the FENE-PG model 
proposed in this study is evaluated by comparison with 
exact Brownian dynamics simulations. 



where, 



(43) 



Since the Rouse matrix is non-singular, the solution to 
Eq. H42() is obtained by setting Uij = 0. This implies, 
^ij", cq ^ ^ij^7 and H* Lj_ cq = 6. With all springs being 
identical at equilibrium , Lj ^q = L^qii, and therefore 



H*Lcq — 1 



(44) 



Using trcr*^ j,q ~ 3, and b* — SA^k. s, it can be shown 



(45) 



using Eq. ^ that for FENE-PG springs, 

3Nl , - ml , + ISA^K, s - 5 



H* = 



3Ni , - ml 3 + 5Ak, s 



For FENE-P springs, Loq = (,q (for all m). where 
was introduced in Eq. H37() earlier. As a result, 



^m, eq 



(46) 



B. FEBS chains with HI 



The choice of parameter values H* and b* for the ap- 
proximate FENE-P and FENE-PG models is made here 
consistently with the procedure used in the case of the 
FENE force law; namely, the equilibrium mean square 
end-to-end distance of a spring predicted by the model 
is equated to bf^N^, s, while the fully stretched length is 
equated to ^kA^k. s- Further, as before. = = 
b'^N^. s/3, and the corresponding time scale As, are used 
as the basic scales for obtaining equations and proper- 
ties in dimensionless form. These assumptions lead to 
^* = Qo/ = 3A^K, s: as in the earlier case. From 
Eq. (|38|l . it is clear that the equilibrium second moments 
a*j predicted by the Gaussian approximation for free- 
draining FEBS chains, are solutions to the following set 



The dimensionless evolution equation for the second- 
moments of a bead-spring chain model with springs obey- 
ing a mean- field potential given by Eq. (|25|l . and with 
fluctuating HI between the beads, is given by, 

H* 

* m=l (47) 

+ A-im • Ljji ■ Qjj^Qj) + —{Aij) , 

Since the tensors X,„ are not direct functions of the 
connector vectors, expectations of the form {Q*Q*a ' Lm ■ 
(C^mj)) appearing in the equation can be simplified by 
assuming that the configurational distribution function is 
Gaussian and using Wick's decomposition rule as in the 
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case of free-draining FEBS chains above. The final equa- 
tion thus obtained for the approximate second-moments 

H* — 

^i], (1) = -r y2 [^trn ■ {Lm ' ^mj + ^mj) 

^ m=l (48) 

+ (A,™ • + A^„J • o-;,^.] + i Ay , 

where, the tensor Aij describes contributions arising 
from averaging the Oseen-Burgers tensor with the non- 
equihbrium distribution function, while the tensor 
represents contributions due to fluctuations in HI. Since 
the expressions for these tensors are lengthy and not very 
illuminating in themselves, they are given by Eqs. (|A1|) 
and (|A2() in Appendix The important point to note 

I 



where, the convention l/VO = has been used. Since 
the modified Rouse matrix is also non-singular, the equi- 
librium solution is once again obtained by setting Uij ~ 
0. Therefore, the equilibrium solution obtained for the 
Gaussian approximation for FENE-PG chains with HI 
is identical to that obtained for free-draining FENE-PG 
chains. In other words, the heuristic two-step procedure 
outlined above to achieve closure in the model combining 
both FE and HI ensures that the equilibrium solution for 
the second-moments is independent of HI. 

Nearly all the closure approximations introduced to 
date for the second- moments ( Q*Q* ) can be derived by 
dropping appropriate terms in Eq. (gHJ. For instance, the 
Zimm model is obtained by (i) using a Hookean spring 
force law, Z„ = 6, (ii) replacing the non-equilibrium 
averaged Oseen tensor contribution Aij with the modi- 
fied Rouse matrix Aij (which arises due to equilibrium 
averaging the Oseen tensor), and (iii) setting = 0. 
Ignoring fluctuations in hydrodynamic interactions in a 
model with Hookean springs ( A^ = and Lm = &) leads 
to the consistent averaging (CA) approximation^, while 
retaining fluctuations in HI leads to the Gaussian approx- 
imation (GA)^. Setting Ay = and L„i = <^,„6 simul- 
taneously (which implies fluctuations in both HI and the 
spring force are ignored), leads to a model combining the 
Consistent- Averaging treatment of HI with the FENE-P 
approximation for the spring force. The predictions of 
this CA-P model for properties in shear flows have been 
studied in detail by Wedgewood and Ottingesi^. On the 
other hand, using only one of the two conditions leads to 
a model in which fluctuations in the corresponding phe- 
nomenon alone arc ignored, leading to the CA-PG model 



here is that Aij and Ay are completely expressible in 
terms of the second moments of the Gaussian distribu- 
tion, and Eq. 1)48(1 is closed with respect to the tensors 

.. . 

At equilibrium, all the tensorial quantities in the equa- 
tions above are isotropic. Furthermore, the HI fluctua- 
tion tensors Ay vanish under isotropic conditions (see 
Appendix^. Consequently, one obtains an equation of 
the same form as Eq. (I42II : 

= ^ [UiraAmj + AimUmj] , (49) 

in which the values of Uij are the same as in Eq. (|43|l for 
the free-draining model with FENE-PG springs, and Aim 
arc elements of the modified Rouse matrix (Eq. ((5()|l ). 



6, (50) 
I 

TABLE I Acronyms used for referring to approximate mod- 
els. 



Acronym 



Treatment of HI 




Free-draining 


FD 


Equilibrium-averaging 


EA 


Consistent-averaging 


CA 


Gaussian approximation 


GA 


Treatment of FE 




Hookean springs 


H 


FENE-P 


P 


FENE-PG 


PG 



or the GA-P model. The model incorporating fluctua- 
tions in both nonlinear phenomena is referred to here as 
the GA-PG model. Note that in all these models, the 
polymer stress is given by Eq. (|39|l . with the appropriate 
choice of Li. 

In discussing predictions obtained using the various 
approximate models, it is convenient to use acronyms 
in which the first two letters refer to the approximation 
used for HI, and the last letters refer to that used for FE. 
Table|l]lists the letter codes used for easy reference. Note 
that in this scheme, the Rouse model is denoted as the 
FD-H model, while the Zimm model is denoted as EA-H. 



V2h* 



9 



C. Diagonalization 



where the auxihary tensors are defined by, 



Previous studies have shown that considerable savings 
in computational time, without a significant loss in accu- 
racy, may be achieved by introducing a diagonalization 
approximationi^'^. The most recent of these approxima- 
tions, carried out in the context of the Gaussian approx- 
imation for Rouse chains with HI, is the "Two-Fold Nor- 
mal" (TFN) approximation of Prakash and Ottinger— 
(named as such because of two different assumptions 
of "normality" -a Normal distribution, and a transfor- 
mation to normal coordinates). Prakash and Ottinger 
showed that for Rouse chains with HI, the predictions of 
the TFN approximations for linear viscoelastic proper- 
ties and steady-state material functions in shear are close 
to those obtained with the original Gaussian approxima- 
tion. Here, we examine the effectiveness of making a 
similar approximation for FENE-P chains with HI. 

The first step in introducing a diagonalization approx- 
imation is mapping the set of the physical connector vec- 
tors {Q* \i = 1, . . . , Ns} onto a set of normal coordinates 
{Q'p\p = Ij ■ • ■ I -^s} using the following transformation. 



Yp, EE H* J2 n.le.A.,H^\ . (58) 



(57) 



Following iPrakash and Ottingeil we assume that the 
normal coordinates covariance matrix is diagonal. 



(59) 



This leads to the following set of equations for the normal 
coordinate variances, 

^'p, K ■ {Xp + Yp) + {Xp + y;) ■ + , 

(60) 



Ns 



where the notation {■■■)p = {■ ■ ■)pp is used. While 
(51) iPrakash and Ottingeil argue that it is not necessary to 
make the additional diagonalization assumptions. 



where H?, are the elements of the orthogonal matrix 
whose columns are formed by the eigenvectors of the 
modified Rouse matrix. 



Ns 



n.ipAyH^g — 6pq Gp , 



(52) 



in which Sp is an eigenvalue of the modified Rouse ma- 
trix. The orthogonal matrix corresponding to the Rouse 
matrix is used for free-draining models without HI, for ef- 
fecting the transformation to normal coordinates. Using 
the orthogonality relation 



Ns 



Ns 



J^^mfe^nfc ~ ^""^ ~ E n|„H|„ , (53) 



fe=l 



^pq ~ ^p^pqi 



XpSpq] 



YpSpq (61) 



in order to derive Eq. H6()|) . Eq. (|55|1 shows that for 
the time derivatives of the off-diagonal components of 
the normal-mode covariances to remain negligible far 
from equilibrium, it is necessary for both Eq. (|59|l and 
Eqs. to be true. 

Equation (|60|l represents a set of ODE's that are 
coupled and nonlinear. The coupling comes about be- 
cause the auxiliary tensors Zp^ Xp and Yp depend on 
all the variances cr*,-. For instance, the FENE-P factors 
appearing in these terms are evaluated as 



1 


1 - 


" Ns 

.P=i 


/b* 



(62) 



the evolution equations for the covariances of the normal 
coordinates ffp^ defined as. 



(54) 



can be shown to be 

N 

i 1 

pg, (1) ^ ^ 4 



CT 



^ Ns 



u=l 

2 ^Pi ' 



(55) 



As a result of the coupling, the CPU-time per time step 
scales as iV^, even though the number of equations is only 
A^s- For every p in Eq. (|60|) . each calculation of Zp, Xp 
and Yp individually involves a double summation over all 
i and j in Eqs. H56|l . H57(l and 1)58(1 . In the calculation of 
the variances using Eq. ((A6() , each of the a*j requires 
a summation over all fc, as can be seen from, 



O", 



Ns 
k=l 



(63) 



Since the values of cr*^ can be stored at the beginning 
of each time-step, the total cost of generating the super- 
matrix containing the a*j tensors scales as . The cost 
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of the calculating S* using 



l—fl 



(64) 



can also be shown to scale as iV^. Thus, as far as the 
CPU-time scaling is concerned, the diagonalization as- 
sumption does not offer an advantage over the original 
approximations in the calculation of properties in un- 
steady flows using bead-spring models with HI. How- 
ever, the reduction in the number of equations to be in- 
tegrated greatly reduces the prefactor in the CPU-time's 
power-law dependence on A^g, and as will be shown in sec- 
tion leads to a significant reduction in CPU time. 

The diagonalized versions of the models will be referred 
to with a "D" prefixed to the acronym representing the 
original approximation, e.g. DFD-P, DCA-P etc. The 
exception to this rule is the TFN-P model, combining 
the TFN approximation for HI with FENE-P springs. 

Before presenting results, the next two sections briefly 
introduce the main material functions examined here, 
and the numerical schemes used to obtain them. 



IV. MATERIAL FUNCTIONS 

Attention is restricted here to some instances of simple 
shear and uniaxial extensional flows. Expressions for the 
material functions in these flows, in terms of components 
of the stress tensor, are discussed below. 



A. Shear flow 

For a simple homogeneous shear flow characterized by 
the shear rate 7, the continuum solvent's velocity compo- 
nents in Cartesian coordinates are given by Vx = 7(i) y, 
Vy = and = 0. The tensor k* = Xg'Vv in Cartesian 
coordinates is 



K*{t* 



^0 1 0> 

7*(r) 1 
,0 0, 



where j*{t*) = 7(<)As is the time dependent dimen- 
sionless shear rate. Here, the time dependence of so- 
lution properties is examined for the sudden imposition 
of steady shear flow: 



nn 




(66) 



While the imposed strain rate 7* at > is a con- 
stant, the solution as a whole is in an unsteady state 
since the properties of the solution change in response to 
the imposition of the flow. Steady-state is reached when 
all properties reach time-independent values which are 



functions solely of the imposed strain rate and the other 
physical parameters characterizing the solution. 

For a simple shear flow defined in Eq. H66|l . the rhe- 
ological properties of interest are the polymer's contri- 
bution to the dimensionless viscosity 77*, the dimension- 
less first normal-stress difference coefficient ^I*^, and the 
dimensionless second normal-stress difference coefficient 
^'2, which are respectively defined by the following rela- 
tions: 



;,* 2 



r — T 



(67) 
(68) 
(69) 



B. Extensional flow 



In a uniaxial extensional fiow with an extension-rate 
£, the velocity field is specified by Vx = i{t)x, Vy = 
-{e{t)/2)y and = -(£(t)/2)z. Therefore, 



'\ 

K*{t*) = e*{t*) \ -i 
,0 



(70) 



where e*[t*) = e{t)Xs is the time dependent dimension- 
less extension-rate. Here, the time dependence of so- 
lution properties is examined for the following kinds of 
unsteady extensional fiows. 

Sudden imposition of steady shear or uniaxial ex- 
tensional flow: 



e*{t*) 



t* <0, 
t* >0. 



(71) 



Sudden imposition followed by cessation of steady 
uniaxial extensional flow: 



(65) 



r < 0, 
o<t* <t* 
t* > ti 



(the stress-growth phase) 
(the stress-relaxation phase) . 

(72) 



The time-dependence of properties in an unsteady state 
is expressed in terms of the strain measure e = e*<*. 

In the case of a uniaxial extensional flow specified by 
Eq. H71|l . the polymer contribution to the dimensionless 
extensional viscosity, denoted as is given by 



— * 
^p 



(73) 



For the fiow pattern described in Eq. (|72|l . rheological 
data is presented in this work in terms of , the magni- 
tude of the polymer's contribution to the dimensionless 
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first normal stress difference: 

NT = It* - T* I (74) 

l,p I p, XX p, yy I V ' / 

We note licrc tliat tfie characterization of a sliear rate 
as "low" or "high" , depends on whether the Weissenberg 
number is much lower, or greater than unity. The Weis- 
senberg number is defined as the product of the shear (or 
extension) rate and a time-scale characterizing the large 
scale dynamics of isolated chains near equilibrium. As is 
well known, this characteristic time-scale is proportional 
to N'^ for free-draining chains, whereas in chains with 
HI, it is proportional to N^^'^ . The Weissenberg number 
in these two cases thus scales as 7 and 7*iV3/2^ 

re- 
spectively. For A'' = 20, a shear rate of 7* ~ 1 therefore 
represents a large strain rate, as will become evident in 
the results below. In extensional flows, the critical strain 
rate for the coil-to-stretch transition (discussed in sec- 
tion marks the boundary between small and large 
extension rates. 



V. NUMERICAL METHOD 

In the case of the un-diagonalizcd approximations, the 
total number of second moments, and therefore the num- 
ber of equations, is equal to 9Ng. However, (T*j = a*/. 
Further, the second moments must not change if the bead 
numbering is reversed. Thus, (7*^ = (y*pj^_i There- 
fore, the number of independent second moments reduces 
by a factor of nearly l/4m. In the diagonalized approx- 
aimations, on the other hand, the number of equations 
for the components of the normal mode variances is just 
9Ns 

The symmetry of the imposed flow field can be used 
to reduce the number of components for all the tensorial 
quantities in the equations. In shear flows, the second 
moment tensors take the following form: 

(CTi (74 \ 
<J5 <J2 . (75) 
(73/ 

In uniaxial extensional flows, in both stress growth and 
relaxation phases, all tensorial quantities in the approxi- 
mations have the form, 

/(Tl 0\ 

a= (72 . (76) 

\ (72/ 

The exploitation of the symmetries above to reduce the 
computation of the fourth rank tensorial function K(s) 
has been discussed thoroughly by Zylka— . In general, for 
any symmetric tensor s, only 21 of the 81 components of 
the K function are non-zero in a transformed coordinate 
system C in which the tensor s has a diagonal represen- 
tation. Out of these 21 components, only 9 require inde- 
pendent numerical evaluation. The components Kim, 



JS'2222 and i^3333 can be used to calculate the remaining 
12 non-zero components. Furthermore, the function K 
always appears in the equations as K: a. In shear flows, 
this means that only 14 components of the K function in 
the laboratory fixed coordinate system are required to be 
evaluated. In uniaxial extensional flows, it can be shown 
that only K^xxx (where x is the direction of stretching) is 
required. As shown by Zylkai^, components of the H{s) 
and K{s) functions in the coordinate system C arc effi- 
ciently evaluated using elliptic integrals. 

The initial condition for the time integration of the 
ode's starting from equilibrium is a*^ = 6ij8 for the 
un-diagonalized approximations. Similarly, in the case of 
the diagonalized approximations, 

"■p.oq = X] n-pCT*,^qH^p = nfp((5y6)n^p = 6 . 

(77) 

As pointed out in the preceding sections, at equilibrium 

-^ij, cq — -^ij^- ^ij, cq — 0; ^ Li.cq ^• 

Since predictions are obtained only for shear and uni- 
axial extensional flows in this study, separate computer 
codes were written for the two different flows, to take 
advantage of the reduction in the components due to 
their individual flow symmetries. A variable-step Adams 
method (the D02CJF routine in the NAG library) is used 
for the integration of the ODE's. 

Steady-state solutions were obtained by setting the 
time derivatives in the ODE's to zero, and then using 
a Wegstein successive substitution algorithm to solve the 
set of nonlinear equations. Solutions are obtained for 
several values of the (shear or extensional) strain rate. 
The converged solution for the set of a*j (or a'^) at any 
strain rate is used as the initial guess for the next higher 
(or lower) strain rate for which a solution is desired. It 
is found that large steps in the strain rates could lead to 
unphysical, non-positive-definite solutions for the second 
momenta^ii^. Such behaviour is avoided by manually 
controlling the size of the steps in strain rate. 

The Brownian dynamics simulations of Rouse and 
FENE chains with (and without) RPY-HI performed 
in this study use a modification of the semi-implicit 
predictor-corrector scheme outlined recently by several 
workeraSi^^ii^, where the error between the predictor and 
corrector steps is used to adaptively control the time-step 
size in every simulation. 

VI. RESULTS AND DISCUSSION 

Results for shear flows and extensional flows are dis- 
cussed separately in the respective sections below. In 
each case, results are presented for (i) models with non- 
linear force laws, but with no HI, (ii) models with HI 
incorporated, but with a linear Hookean force law, and 
(iii) models with non-linear force laws and HI incorpo- 
rated. In particular, in case (i), results obtained with 
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the new FENE-PG force law introduced here are com- 
pared with results obtained with the FENE-P force law 
and with BDS, in case (ii) (which involves only models 
that have been introduced previously) , results for certain 
material functions that have not been predicted earlier, 
are presented, and in case (iii), the combined influence 
of both phenomena are examined, aided by the results 
obtained in cases (i) and (ii). 



A. Shear flow 

1. Startup of shear flow 

The performance of the FENE-P and FENE-PG ap- 
proximations for free-draining FEBS chains in the pre- 
diction of unsteady behaviour in shear flows is first con- 
sidered. At low shear rates close to equilibrium, finite 
chain extensibility has little influence on the dynamic be- 
haviour predicted by the models. Therefore, the predic- 
tions of the free-draining FEBS models are close to those 
predicted with the Rouse model. At higher shear rates, 
earlier studies with free-draining FENE-P dumbbells^ 
and chains with N = 10^ show that this approximation 
leads to predictions for unsteady shear properties that 
agree only qualitatively with the exact results obtained 
using FENE springs. Figures ^ (a-) and (b) show that 
this trend is continued for 20-bead chains. The incorpo- 
ration of fluctuations in the spring force with the use of 
the FENE-PG approximation leads to a clear improve- 
ment in the quality of the predictions for all properties 
at moderate shear rates (represented by the curves for 
7* = 0.23). As the dimensionless shear rate approaches 
unity, however, the FENE-PG approximation also begins 
to deviate from the exact BDS results. Interestingly, it is 
seen that the FENE-PG force law leads to the prediction 
of a transient non-zero positive second-normal stress dif- 
ference coefficient which approaches zero at steady-state, 
in good agreement (within errors in the simulations' re- 
sults) with the results of simulations. 

A striking feature is the large overshoots predicted by 
all the models in the growth of rj* in Fig. (a). A 
similar prediction is obtained for This feature is 

not observed at any shear rate for Rouse chains with HI 
(see Fig. |21 (a)), and neither is it observed at low shear 
rates for FEBS chains. The overshoots in the material 
functions appear to be accompanied by a corresponding 
overshoot in the mean-squared end-to-end distance of the 
coils. Overshoots in the transient i]*, and i?*^ during 
start-up of shear flows are kn own to be a signature of 
the influence of iHur et aP have obtained 

experimental evidence of such overshoots using dilute 
(and semi-dilute) solutions of DNA, and have shown us- 
ing simulations of bead-rod and bead-spring chains that 
the overshoot is caused by the rapid localized stretching 
of the molecules by the flow, before the orientation of the 
chains towards the flow direction reduces the total solvent 
drag sufficiently to reduce the stretching. Herrchen an 
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^K,s = 18-3 
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FD-PG 



Y = 0.91 




FIG. 1 Growth of (a) viscosity and (b) second normal-stress 
difference coefficient, for free-draining FEBS chains during 
start-up of steady shear fiow. 



dOttinger— observed that the overshoot predicted with 
FENE-P springs is larger than that obtained with FENE 
springs in the case of dumbbells. 

Although the qualitative behaviour of the different ap- 
proximations for HI in Rouse chains has been well docu- 
mented in the case of steady shear flows, their predictions 
of rheological and conformational properties in transient 
shear flows have not received much attention. For Rouse 
chains with HI, use of equilibrium averaging in the Zimm 
models results in a fixed time dependence of properties 
that is independent of the shear ratai^. This behaviour 
parallels the shear rate independence of the material 
functions predicted at steady-state by the Zimm model. 
In contrast, the predictions of the consistent-averaging 
and Gaussian approximations, and the exact results of 
the simulations for the growth of the shear material func- 
tions depend on the imposed shear rate, as demonstrated 
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FIG. 2 Growth of (a) viscosity, and (b) second normal-stress piQ 3 Growth of (a) first normal-stress difference coefficient, 
difference coefficient, for Rouse chains with HI during start-up and (b) viscosity, for FEES chains with HI during start-up of 
of steady shear flow. steady shear flow. 



in Figs. 121(a) and (b) for the transient variation in 77* and 
^'2, respectively, after the imposition of a steady shear 
rate. 

In Fig. 121(a), the curves predicted with the Gaussian 
approximation initially lie below those predicted with 
consistently averaged HI, for both the shear rates con- 
sidered. At larger values of t* , as steady-state is ap- 
proached, the prediction of the Gaussian approximation 
at the higher shear rate is seen to overtake the CA-H 
curve, whereas at the lower shear rate, the GA-H predic- 
tion is always below that predicted by the CA-H model. 
The importance of including a description of the fluctua- 
tions in approximate models accounting for HI is clearly 
highlighted by the superiority of the predictions of the 
Gaussian approximation. Interestingly, Fig. (21 (b) shows 
that fluctuations in HI lead to a negative second nor- 
mal stress difference coefficient soon after the inception of 
flow. Zylka— had earlier demonstrated the accurate pre- 



diction of the sign of with the Gaussian approximation 
at steady state. The results here show further that the 
approximation is excellent even in transient shear flows. 

Figure 121(a) shows that the GA-PG approximation is 
particularly useful in shear flows of moderate strength 
where fluctuations in the spring forces and HI are both 
important. The prediction of the GA-PG model is in re- 
markable agreement with simulations' results for the un- 
steady growth in the first normal-stress difference coeffi- 
cient following the imposition of a moderate shear rate of 
7* = 0.23. Similar quantitative agreement is obtained for 
the unsteady variation in 77* and i?* ^ at this shear rate. 
At lower shear rates, FE is not important, and the be- 
haviour is the same as that obtained with Rouse chains, 
for which it was shown above that the Gaussian approx- 
imation does well in describing unsteady behaviour. 

At high shear rates, where FE becomes the dominant 
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FIG. 4 Comparison of the predictions of the GA-H, FD-PG 
and GA-PG models for the growth of the shear viscosity with 
BDS results obtained with Rouse chains with HI (green trian- 
gles), free-draining FENE chains (grey squares), and FENE 
chains with HI (magenta circles). The dimensionless shear 
rate is unity for Rouse chains and 0.91 for FENE chains. 



effect, however, predictions of no single model are uni- 
formly in agreement with the BDS results for all the 
properties examined, as illustrated in Fig. 13 (b). To un- 
derstand the source of the deviations in the predictions 
of the GA-PG model, we compare in Fig. 0] predictions 
of the GA-H, FD-PG and GA-PG models at 7* ~ 1. 
The transient evolution of the viscosity in the GA-PG 
model can be broken into two phases: the early over- 
shoot phase, and the later steady-state approach phase. 
At early times, we see that the GA-H model is in very 
good agreement with the BDS data (green triangles) for 
Rouse chains with HI. In the same period, the FD-PG 
predictions deviate quite strongly from the corresponding 
BDS data for free-draining FEBS chains (grey squares). 
The agreement of the composite GA-PG model with the 
BDS data for FEBS chains with HI (magenta circles) is 
thus seen to be intermediate in the overshoot phase. To 
understand the behaviour at steady state, we note that 
fluctuations in HI at high shear rates aid in the pene- 
tration of the solvent's velocity field into the deformed 
polymer coil, and the viscosity predicted is higher and 
closer to predictions with the corresponding free-draining 
model. This important aspect of Hi's influence on the be- 
haviour of dilute polymer solutions will be more clearly 
demonstrated shortly in predictions of steady state prop- 
erties. From the overprediction of the viscosity by the 
GA-H model at large times in Fig. 0] we see that at high 
shear rates, the Gaussian approximation for HI overesti- 
mates the influence of fluctuations in HI. In Fig.|21(b), we 
see that the CA-PG model's prediction is lower than the 
BDS results. The Gaussian approximation for HI thus 
pushes the viscosity prediction of the CA-PG past the 
BDS results. 
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FIG. 5 Variation of (a) the steady-state first normal-stress 
difference coefficient, and (b) the steady-state mean-squared 
end-to-end distance, with shear rate, for free-draining FEBS 
chains. 



2. Steady shear flow 

The steady-state shear viscosity predicted by the 
FENE-P approximation has been observed earlier to be 
in quantitative agreement with the values obtained with 
BD simulations of small free-draining FENE chains^LS. 
For chains with N ~ 20, it is observed that the pre- 
dictions of 77* obtained with the FENE-P and FENE- 
PG approximations are nearly identical. However, in 
Figs. El (a) and (b), it is seen that the FENE-PG ap- 
proximation for the FENE force law does much better 
than the FENE-P approximation in its prediction of the 
steady-state Vl** and Rl'^, respectively, at moderate shear 
rates. At 7* « 1, the predictions of i?*^ and '^l with 
FENE-PG springs also begin to deviate from BDS results 
significantly, but are closer to the exact results than the 
predictions obtained with FENE-P springs. From the 
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fact that the stresses and R* ^ values predicted with sim- 
ulations and the FD-PG model are lesser than the predic- 
tions of the FD-P model; we conclude that fluctuations 
in the spring forces lead to an increase in the effective 
stiffness of the chains. 

As mentioned before, the predictions of the Gaussian 
approximation for Rouse chains in steady shear flows 
have been compared earlier with those of the Consistent- 
and Equilibrium- Averaging (Zimm) approximations, and 
with the results of exact BD simulations^. Figure El il- 
lustrates some of the key features of the predictions of 
the different approximations. As is well known, the net 
reduction in the average solvent velocity gradient within 
the polymeric coili22il2i^ caused by HI results in smaller 
values of the shear viscosity at low to moderate shear 
rates, when compared to the free-draining Rouse model. 
The fact that the predictions of the Gaussian approx- 
imation and BD simulations are lower in Fig. than 
those obtained with the Zimm and consistent-averaging 
approximations near equilibrium and at moderate shear 
rates, indicates that under these conditions, fluctuations 
in configurations at equilibrium result in a more effec- 
tive screening of the velocity field within the coil. As 
the shear rate increases, the steady-state viscosity and 
the first normal stress difference coefficient first decrease 
in models with HI, other than the Zimm model. This 
shear-thinning, although not very strong, is a distinc- 
tive feature of configuration-dependent HI. The initial 
shear-thinning is observed to be more prominent with 
consistently-averaging than with the Gaussian approxi- 
mation for HI. As the shear rate increases, the predic- 
tions of 77* and '^l obtained with BD simulations, and 
the CA-H and GA-H models, reach a minimum with re- 
spect to the shear rate 7*. With a further increase in the 
shear rate, shear-thickening is observed in these models, 
wherein 77* and begin to increase towards the con- 
stant value predicted by the Rouse model in which HI 
is absent. The increase in 77* and towards the Rouse 
model's predictions is understood to be the consequence 
of the decreasing influence of HI as the average separa- 
tion between the beads in the chains increases. Interest- 
ingly, in Fig. El the predictions of the GA-H model and 
the BD simulations' results lie above the curve predicted 
by the CA-H model beyond a threshold shear rate, and 
are closer to the free-draining Rouse model's prediction. 
This behaviour is contrary to that observed near equi- 
librium, and shows that the role played by fluctuations 
in HI reverses at higher shear rates, and in fact aid the 
penetration of the solvent velocity field into the polymer 
coil as the coil deforms and becomes highly anisotropic. 

It is known that both the consistent-averaging and 
Gaussian approximations lead to shear rate dependent 
predictions for ^Pj, which are opposite in sign for a large 
range of shear rates. This was first shown by Zylka^, thus 
demonstrating that fiuctuations in HI lead to a change in 
the sign of at low to moderate shear rates for Rouse 
chains, inline with the results obtained here above for 
unsteady shear fiow. 
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Figure [71 shows the results for the steady-state ^* ob- 
tained with the different approximate models combining 
HI and FE. Predictions for the transient growth of 
at 7* = 0.23 were presented earlier in Fig. (a). The 
qualitative features of the predictions for the steady-state 
viscosity r/* are the same, although the shear-thinning- 
thickening behaviour is more pronounced in the case of 
^'*. As before, the predictions of these models are com- 
pared with the results of BD simulations of FENE chains 
with fluctuating HI. 

It is clear from Fig. [3 that among all the approximate 
models, only the predictions of the GA-PG model are 
close to the simulations' results across the whole range 
of shear rates examined. As one would expect, at low 
shear rates where the influence of FE is negligible, the 
observed behaviour is close to that observed for Rouse 
chains with HI, for which the Gaussian approximation 
for HI has been shown earlier to perform better than 
the Consistent-Averaging approximation. At moderate 
shear rates when FE begins to exert its influence, it is 
seen that the Gaussian approximation is not as accurate 
when used with FENE-P springs as it is with FENE- 
PG springs. This suggests that fluctuations in the spring 
forces also become important at moderately large shear 
rates. Interestingly, spring force fluctuations are how- 
ever observed in Fig. |S1 to have a negligible effect on the 
steady-state ^fj, and the predictions with the FENE-PG 
force expression are nearly the identical to those obtained 
with the FENE-P approximation. 

As the shear rate increases further, the differences be- 
tween the consistent-averaging and Gaussian approxi- 
mations reduce for the same choice of the spring force 
law. However, it is observedSi that at high shear rates, 
the predictions with models using the Gaussian Approx- 
imation for HI are closer to their free-draining counter- 
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FIG. 7 Variation of the steady-state first normal-stress dif- 
ference coefficient with shear rate, for FEBS chains with HI. 
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FIG. 9 Prediction of shear-thinning-thickening-thinning of 
the steady-state shear viscosity by models with HI, for TVk = 



Approximation for HI and the FENE-PG approximation 
for the nonlinear spring force is approximately bounded 
from above by the envelope formed by the Gaussian Ap- 
proximation for Rouse chains at low to moderate shear 
rates, and by the terminal shear-thinning curve predicted 
by the free-draining model with the FENE-PG springs at 
high shear rates. The shear rate at which the 77*-versus- 
7* curves predicted by the two simpler models (viz. the 
GA-H and FD-PG models) intersect can be regarded as 
a critical shear rate 7*. At shear rates smaller than 7*, 
the overall behaviour of the polymer chains is dominated 
by HI, whereas EE is the dominant phenomenon when 

7* > 7c- 
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FIG. 8 Variation of the steady-state second normal-stress 
difference coefficient with shear rate, for FEBS chains with 
HI. 



B. Extensional flow 

1. Start-up and cessation of steady extensional flow 



parts, than those using consistently-averaged HI. This 
behaviour once again points to the fact that fluctuations 
in HI aid in the penetration of the solvent's velocity gradi- 
ent into the deformed polymer coil at higher shear rates. 

Although at high shear rates Rouse chains with HI 
predict shear-thickening, FEBS models with HI predict 
shear-thinning. Kisbaugh and McHughi^ showed using 
the diagonalized version of the CA-P model that com- 
plicated shear-tinning-thickening-thinning behaviour can 
occur as a result of the interplay between the emergence 
of EE and the fading of HI. Figure El shows the predic- 
tions of ry* obtained with the CA-PG and GA-PG models 
for chains with 25 springs, and TYk, s = 400. The be- 
haviour in Fig. |51 suggests that the curve for the steady- 
state ?7* predicted by a model combining the Gaussian 



Models using Hookean springs predict unbounded ex- 
tension is strong extensional flows. As a result of this 
unrealistic behaviour, the predictions of approximations 
for HI in Rouse chains have not been examined in strong 
and unsteady extensional flows. However, at any finite 
time following the imposition of the flow, the values of the 
stresses and other properties obtained with Rouse chains 
are still finite. Therefore, by comparing the predictions 
of the approximations and the results of BD simulations 
for Rouse chains one can study in isolation the influence 
of HI in strong extensional flows. 

Figure ^| shows the growth of TVj* p upon the imposi- 
tion of a steady extensional flow of e* = 0.04 (which is 
greater than the critical strain rate for the coil-stretch 
transition el = 0.012 for iV = 20 and h* = 0.25). Also 
shown in Fig. ^| are the predictions for the relaxation of 
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FIG. 10 Growth and decay of the polymer contribution to 
the first normal-stress difference for Rouse chains with HI, 
during start-up and following cessation of steady extensional 
flow, respectively. 
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FIG. 11 Effect of fluctuations in HI on the growth and decay 
of the polymer contribution to the first normal-stress differ- 
ence for Rouse chains with HI, during start-up and following 
cessation of steady extensional flow, respectively. 



N* p after the cessation of the extensional flow at differ- 
ent values of the Hencky strain. 

The results in Fig.llOlshow several fascinating features. 
Firstly, the predictions of all the approximate models, 
including the Zimm model with equilibrium-averaged HI 
appear to be in remarkably close agreement with the ex- 
act results of BDS in the stress growth phase. The Zimm 
model only begins to deviate noticeably from the rest of 
the data during stress relaxation after flow stoppage at 
high strains. Similar results are also obtained for i?*^— . 
Although the differences between the different approxi- 
mations appear to be small on the scale of the plot in 
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FIG. 12 Growth of (a) the polymer contribution to the first 
normal-stress difference, and (b) the mean-squared end-to- 
end distance of free-draining FEES chains, during start-up, 
and their relaxation following cessation, of steady extensional 
flow. The horizontal red line in (b) indicates the maximum 
possible value of i?* = SA'k. sN^. 



Fig. El in Figure im the differences between the predic- 
tions of these approximations are examined more closely 
by plotting the relative difference between the predictions 

of Nip of the two approximations. 
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It is seen that fluctuations in HI can account for as much 
as 30% of the polymer stress in extensional flows. 

Predictions of free-draining FEES models employing 
the FENE-P and FENE-PG approximations are com- 
pared with exact results of the simulations for the growth 
and relaxation of ^ and R*'^ in Figs. El (a) and (b), 
respectively. Also shown for comparison are the predic- 
tions of the simple Rouse model. The strain-rate for 
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which the data are shown in Fig. ^] is larger than the 
critical strain rate i* = sin^(7r/2A^) = 6.2 x 10~^ for 
free-draining chains with — 20. At this relatively high 
strain rate, the Rouse model predicts an exponential in- 
crease in both Nip and R*^ during extension. At small 
values of the Hencky strain during start-up of the exten- 
sional flow at this strain rate, all the data for the dif- 
ferent models are essentially identical, since the springs 
are not stretched sufficiently for the spring-force nonlin- 
earities to have any signiflcant effects. The influence of 
FE is first discerned in the growth of the stress when the 
curves for the FEBS chains separate from the prediction 
of the Rouse model. At higher strains, the deviations be- 
tween the approximations and the exact results for A^^* p 
become more marked and reach a maximum when the ap- 
proximate results level off sharply towards their eventual 
steady-states. The approach of the BD simulations' data 
to steady-state is more gradual. In the stress growth 
phase, the predictions of the FENE-PG approximation 
are closer to the exact results, than those obtained with 
the FENE-P approximation. 

The relaxation behaviour of free-draining FEBS chains 
following the cessation of extensional f low has been previ- 
ously studied^i^. While iDovle et al.l compared the pre- 
dictions of the FENE-PM approximation with BD sim- 
ulations of FENE dumbbells and also bead-rod chains, 
IWiesti examined some interesting qualitative features in 
the stress-relaxation phase using the FENE-P model. 

The FENE-PG approximation is seen to lead to more 
accurate description of both the stress and extension in 
the relaxation phase when the strain at flow stoppage is 
not much larger than the strain at which FE begins to 
exerts its influence, whereas the FENE-P approximation 
leads to an overprediction of both properties. When the 
Hencky strain at flow stoppage is larger, the FENE-PG 
force law is seen to be in fact less trustworthy. This loss in 
accuracy of the FENE-PG model when chains are highly 
stretched is consistent with the behaviour observed ear- 
lier in steady and unsteady shear flow. 

The predictions of approximate models with both non- 
linear phenomena incorporated, for the unsteady growth 
and relaxation of A^^*p and R*^ following start-up and 
cessation of an extensional flow are shown in Figs. ll3l fa') 
and (b). The initial behaviour {e < 4) in the stress- 
growth phase is close to that observed for Rouse chains 
with HI since the influence of FE is weak at these strains. 
The differences between the various approximations grow 
larger at intermediate strains during extension. At these 
strains and above, the curves predicted by the combined 
models are qualitatively similar to the predictions ob- 
tained with the approximations for free-draining FEBS 
chains. In particular, the nature of deviations between 
the BD simulations' results and the approximations' pre- 
dictions appear to be decided largely by the approxima- 
tion used for the FENE nonlinearity. At higher strains, 
as chain become highly stretched, incorporating fluctua- 
tions in HI through the Gaussian Approximation make 
the curves predicted by the combined models lie closer 
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FIG. 13 Growth and decay of the polymer contribution to (a) 
the dimensionless first normal-stress difference, and (b) the 
mean-squared end-to-end distance, for FEBS chains with HI, 
during start-up and following cessation, of steady extensional 
flow. The horizontal red line indicates the maximum possible 
value oi Rl^ = iN^,sNi. 



to those obtained with the corresponding free-draining 
modelSi. Predictions of the relaxation behaviour are ob- 
served to be better with the FENE-PG force law when 
the strain at flow stoppage is low. The largest deviations 
in the approximations' predictions of the relaxation be- 
haviour occur when flow is stopped at a strain in the 
"knee" region just prior to steady-state. On the whole, 
the GA-PG model appears to be quite accurate in its de- 
scription of the normal stresses and chain deformation in 
unsteady extensional flows. 

The prediction of stress-conformation hysteresis is 
another test of the accuracy of closure approxima- 
tions. Stress-conformation hysteresis in dilute poly- 
mer solutions has been studied in experimentsSiS and 
theoreticalljiSi^SiiSiSii^. In extensional flow experi- 
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FIG. 14 Stress-conformation hysteresis for FEES chains with 
and without HI. 



merits, it is observed that when the normal stress is plot- 
ted against the birefringence of the solution, the curves 
described by the data points during the stress-growth 
and stress- relaxation phases are very different. In par- 
ticular, at any fixed value of the birefringence, the stress 
during extension is larger than that measured during re- 
laxation. Conversely, at a fixed stress, the birefringence 
during relaxation is larger. 

Figure ^1 compares the A^j* p- versus- An* behaviour 
predicted by BD simulations and the various approxi- 
mations. It is immediately clear that closure approxi- 
mations fail to reproduce the size of the hysteresis loops 
observed in the simulations. Models using the FENE-P 
springs predict larger birefringence than the exact simu- 
lations at any level of the stress. When fluctuations are 
accounted for in the spring forces through the FENE-PG 
force law, the size of the hysteresis loop becomes even 
smaller. Nevertheless, the curves predicted with FENE- 
PG springs now lie within the loops predicted by the BD 
simulations. 

Remarkably, HI appears to have almost no influence 
on stress-configuration hysteresis. Although the maxi- 
mum values of the stress and An* are larger when free- 
draining models are used, the simulations' and approx- 
imations' predictions with and without HI show signifl- 
cant overlap. Similar behaviour is also observed on a plot 
of A''* p-versus-i?* ^— . The reason behind this insensitiv- 
ity of stress-conformational hysteretic behaviour to HI is 
not understood. 



2. Steady extensional flow 

Figure [Tsl compares the predictions of the steady-state 
polymer extensional viscosity by different approxima- 
tions with the simulations' results obtained for FEBS 
chains with and without HI. The BDS results were ob- 
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FIG. 15 Variation of the steady-state polymer extensional 
viscosity with extension-rate, for FEBS chains with HI with 
and without HI. The red square and circle symbols represent 
results obtained with BDS of FEBS chains with and without 
HI respectively. The horizontal dashed red line is the dimen- 
sionless extensional viscosity of fully-stretched free-draining 
chains calculated using Eq. 1791 and marks the upper limit 
of the extensional viscosity for FEBS chains with parameter 
values shown. 



tained using step changes in the extension rates and in- 
tegrating until a steady value is obtained for 77* at each 
extension rate. The most important characteristic of 
the curves for the steady-state extensional viscosity is 
the occurrence of a sudden rapid increase in its value 
within a narrow range of strain rates around a criti- 
cal elongation rate e*. As is well known, the value of 
e* can be calculated analytically in the Rouse model to 
be e* = 0.5/ti = sin2(7r/2iV)i& at which the steady- 
state 77* predicted by the Rouse model diverges. A sim- 
ilar divergence is predicted with the Zimm model at 
e* = O.S/rj'^. Here, rj*"' and tI'^ are the largest time 
constants is the respective spectra of relaxation times in 
the Rouse and Zimm models. For FEBS chains with or 
without HI, 77* remains bounded by the limit for fully- 
stretched free-draining chains completely aligned in the 
direction of extension^, 
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As observed earlier in shear flows, the steady-state 77*- 
versus-e* curves predicted by models with FEBS chains 
are indistinguishable from those obtained with Rouse 
chains until FE becomes important. Predictions with 
the FD-P approximation in the super-critical regime 
(e* > £*) are in remarkable agreement with the BDS 
results, whereas the FD-PG approximation significantly 
underpredicts the exact results. The accuracy of the FD- 
P approximation in its prediction of the steady-state ex- 
tensional viscosity is well knownS^. This behaviour how- 
ever is in contrast with that observed in steady shear 
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flows (Fig.O, where the FD-PG approximation's predic- 
tion are closer to the simulations' results even at high 
shear rates. 

As mentioned earlier, Lielens et alm^ and Wies^ have 
argued that the accuracy of the FENE-P approximation 
at high extension rates stems from the fact that the exact 
distribution of spring lengths is very narrow and close to 
being a (5-function centered at some value Q lesser than 
Qq. Under such conditions, the FENE nonlinearity ^ is 
well approximated by the FENE-P function ^. When Q 
is not very close to Qq, it is clear that the additional 
term in the FENE-PG function L in Eq. (EUl leads to 
inaccuracies. 

We recall that higher order terms in the series expan- 
sion in Eq. (|28|l have been dropped in obtaining L in 
Eq. H34|) . Since every term [Eq. (|30|l ] in the series in 
Eq. (|28|l is positive, it is evident that after applying the 
Peterlin approximation to every term in the series, the 
full summation (if convergent) will lead to an even greater 
overestimation of the influence of spring-force fluctua- 
tions when the true spring length distribution is narrow. 
Since the series expansion in Eq. H28|) itself has been de- 
rived under the assumption of a Gaussian distribution, it 
appears that the extension of the Gaussian approxima- 
tion to the FENE springs may not be appropriate in situ- 
ations where the spring length distribution is narrow, but 
Q is not close to Qq. The better agreement of predictions 
with FENE-PG springs with BDS data observed earlier 
in shear flows not only indicates that spring length fluc- 
tuations are larger in shear flows, but also shows that the 
Gaussian approximation captures the influence of these 
fluctuations more accurately. As the strain-rate increases 
and Q approaches Qo, it was shown in section FlILAI that 
the additional contribution in the FENE-PG spring force 
becomes negligible in comparison with ^, and the differ- 
ences between the predictions of the two approximations 
become smaller. 

For chains with HI. the superiority of the Gaussian 
approximation is clearly evident in Fig. 1151 in the sub- 
critical regime. At strain-rates greater than e*, it is 
the GA-P model that performs best (also sec Fig. [TBI 
below). This is somewhat surprising, since one might 
expect from the better performance of the FD-P model 
for free-draining chains observed earlier, that fluctuations 
are small when chains are stretched out in the super- 
critical regime. Therefore, if the influence of both spring- 
force and HI fluctuations are unimportant, the predic- 
tions of the CA-P model should be close to the BDS data. 
Instead, however, the GA-P model does much better than 
the CA-P model in Fig. ^1 One plausible explanation 
for this observation could be that although fluctuations 
in spring lengths may be small and the FENE-P approxi- 
mation for the spring force is accurate, the distribution of 
overall chain configurations may not be as sharply peaked 
at strain rates above e*. The resultant fluctuations in 
HI may therefore may still play an important role. The 
GA-PG model underpredicts the BDS results for the ex- 
tensional viscosity when e* > i^, and this is seen to be 
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FIG. 16 Prediction of coil-stretch hysteresis in the steady- 
state extensional viscosity by closure approximations for 
FEBS chains with HI. The upright and inverted triangle sym- 
bols represent BDS data obtained by successively stepping up, 
and stepping down, the extension rate, respectively. 



the result of the inaccuracy of the FENE-PG spring force 
at high extensions. 

It is observed that the approximations for FEBS chains 
with HI can predict multiple steady-state values for the 
extensional viscosity for a range of values of the exten- 
sion rate e* [Fig. ^]. The variation of the steady-state 
extensional viscosity with extension-rate predicted by the 
approximations in Fig. llBl was obtained as follows. Start- 
ing with the equilibrium second moments (cr*j = SijS) 
as the initial guess, the set of coupled, nonlinear equa- 
tions for the steady-state second moments were first it- 
er atively solved at a small value (<C e*) of the extension 
rate. After calculating the extensional viscosity with the 
converged solution for the equilibrium second moments, 
this solution is used as the guess for obtaining the so- 
lution numerically at a slightly larger value of e*. This 
process is continued until e* is well beyond the value of e* 
for the coil-to-stretch transition. The ry*-versus-£* curve 
obtained thus is indicated as having been obtained with 
a "coiled initial condition" in Fig.^1 The process is then 
reversed by successively lowering the extension-rate. In 
Fig. El it is seen that with the notable exception of the 
EA-P model; the predicted r7*-versus-£* curve follows a 
different route as e* is decreased for all the other approx- 
imations for chains with HI. Since this steady-state solu- 
tion branch is obtained starting from a highly "stretched 
initial condition" , it is indicated as such in Fig. E| A 
similar procedure described earlier by Schroeder et alm^ 
was used to obtain the BDS data. 

The existence of such "coil-stretch hysteresis" was first 
proposed in a landmark study by De GennesS, who used 
a variable friction coefficient in a simple dumbbell model 
to model the effect of the deformation of the polymer 
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coil on its hydrodynamics. A further closure approxima- 
tion equivalent to the consistent-averaging treatment was 
used to obtain a solvable dumbbell model that predicted 
coil-stretch hysteresis in a uniaxial extensional flow. To 
the best of our knowledge, Fig.^]is the first demonstra- 
tion that closure approximations for FEBS chain models 
with HI can predict coil-strctch hysteresis. 

It is interesting to note that the size of the window 
of strain rates where multiple steady-states are obtained, 
depends strongly on the approximation used. The win- 
dow size appears to depend on the difference between 
the extensional viscosities in the coiled state (the lower 
branch of the hysteresis curve) and the stretched state 
(the upper branch of the hysteresis curve). Fluctuations 
in HI lead to a widening of the window because these 
fluctuations decrease ^* in the coiled state, but increase 
77* when the chains are highly stretched. As before, the 
BDS data are quite close to the predictions of the GA- 
P model. Significantly, multiple steady-states are never 
obtained with the EA-P m odel, and in the free-draining 
models. As pointed out by iDe Genne^ . coil-stretch hys- 
teresis occurs principally because of the modification of 
the mobility matrix due to the deformation of the poly- 
mer coil. With equilibrium averaged HI (and in all free- 
draining models), the mobility matrix is deformation- 
independent, and hence no coil-stretch hysteresis is ob- 
served. 

As mentioned above, De Gennes first obtained such 
multiple steady-states using a mean-field closure approx- 
imation in a simple dumbbell model with a variable 
configuration-dependent drag coefficient. Later, Fan et 
alm^ obtained the solution to the Fokker-Planck equa- 
tion for the dumbbell model numerically without using 
any closure approximations, and found that the resulting 
steady-state expectations are single-valued functions of 
the extension rate. They argued that De Gennes' predic- 
tion of a coil-stretch hysteresis must therefore be an arti- 
fice of the mean-field approximations used in his model. 
More recently, Schroedcr et alm^ showed that the exact 
steady-state configurational distribution in the variable 
drag dumbbell model is bi-modal, and that a large (effec- 
tive) free energy barrier between the coiled and stretched 
states leads to ergodicity-breaking. In contrast, however, 
the use of closure approximations leads to two distinct 
uni-modal Gaussian solutions to the modified Fokker- 
Planck equations. 

The BDS results obtained by Schroeder et alm^, and 
those in Fig.lTCIsuggest that the exact distribution of the 
steady-state extensional stress in FEBS chains with HI is 
also bi-modal. The linearity of the Fokker-Planck equa- 
tion [Eq. 1^] in the configurational probability distribu- 
tion function further guarantees that its exact solution is 
unique. The introduction of any mean-field approxima- 
tion however leads to a modified Fokker-Planck equation 
that is nonlinear in the probability distribution function, 
since the coefficients in the modified Fokker-Planck equa- 
tion are functionals of the distribution function. While 
the solution of this modified equation may be a uni-modal 



Gaussian distribution, the nonlinearity in the distribu- 
tion function means that the solution may no longer be 
unique. Thus, the exchange of the nonlinearities due to 
FENE and HI for a nonlinearity in the distribution func- 
tion destroys the multi-modal nature of the original so- 
lution, but allows the modified Fokker-Planck equation 
to have multiple solutions. 

In spite of the fundamentally different probability dis- 
tributions, the results presented in Fig. ^| above indi- 
cate that the multiple steady-states obtained with the 
approximations can closely follow the long-lived kinet- 
ically trapped states in the original model. Therefore, 
closure approximations for FEBS chains with HI may 
prove to be useful in a detailed exploration of hysteretic 
phenomena caused by ergodicity-breaking in dilute poly- 
mer solutions. 

We next briefly discuss the accuracy of the diagonal- 
ization approximations in shear and extensional flows. 



C. Two-fold Normal Approximation 

Figure [T7I compares the predictions at high shear and 
extension rates obtained with the diagonalizcd DGA-H, 
DFD-P and TFN-P approximations and their respective 
un-diagonalized versions, namely, the GA-H, FD-P and 
GA-P approximations. The effect of making the diag- 
onalization assumption is examined here by calculating 
the relative deviation between the predictions of a prop- 
erty TT obtained with an approximation "A" and its di- 
agonalizcd version "DA" : 

A,,.^^^H£^. (80) 
tta 

Figures El {^) a-nd (b) demonstrate that the relative de- 
viation between the predictions of GA-P model and its 
diagonalizcd version in shear flow, the TFN-P model, can 
be quite significant. In contrast, the diagonalization as- 
sumption is much more accurate for approximations for 
Rouse chains with HI (the GA-H and TFN-H models), 
for which the largest relative deviation does not exceed 
10% in Fig. ini (a) and (b). On the other hand, the 
overall time variation of the deviation between the TFN- 
P and GA-P models is qualitatively similar to that ob- 
served between the diagonalizcd FD-P (DFD-P) model 
and the FD-P model, but is smaller in magnitude. The 
magnitude of the deviations depend on the property con- 
sidered. In shear flows, the deviations in the predictions 
for rj* are smaller. In both shear and extensional flows, 
the deviations are larger in the predictions of R*^ than 
in the stresses. 

The results above indicate that for Rouse chains, or 
for FEBS chains at low strains and/ or strain-rate, when 
EE is unimportant, the inclusion of HI diminishes the 
influence of the off-diagonal components of the normal- 
mode covariances on the overall behaviour predicted by 
the models, and the diagonalization approximations are 
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FIG. 17 Predictions of the diagonalization approximations 
for (a) the growth of the first normal-stress difference coef- 
ficient after imposition of a steady shear rate, and (b) the 
growth and relaxation of the polymer contribution to the 
dimensionless first normal-stress difference, during start-up 
and following cessation, of steady extensional fiow. In (a), 
iVK, s = 18.3 for models with FENE-P springs, and h* = 0.25 
for models with HI. In (b), Nk, s = 50 for models with FENE- 
P springs, and h* = 0.25 for models with HI. The insets 
plot the relative deviations of the predictions of the diago- 
nalized and un-diagonalized approximations, calculated using 
Eq. JHUJ. 



quite accurate. The diagonalization assumption begins 
to break down when FE begins to exert its influence. 

In spite of the deterioration in the quality of the di- 
agonalization approximations when FE is important, the 
massive reduction in CPU-time usage (Table achieved 
with these approximations make them an ideal tool for 
rapid preliminary exploration of vast regions of parame- 
ter space. In fact, Eq. (|6()() coupled with the equation for 
the polymer stress in Eq. H39(l , constitute a new constitu- 
tive model for the polymer stress in dilute polymer solu- 



TABLE II Sample comparison of CPU-time requirements 
for BD simulations, un-diagonalized and diagonalized approx- 
imations. The value reported for BD simulations is for 10'* 
trajectories. For A'^ = 160, the reported value has been calcu- 
lated from the CPU-time required for 200 trajectories. The 
CPU-times reported for the CA-P and GA-P approximations 
for N = 160 have been projected from the CPU-time-versus- 
N scaling observed for 20 < A < 100. All calculations were 
performed on 2.66 GHz, Pentium 4 (Dell Precision 250) pro- 
cessors. 



A 




CPU-time in hrs. 






BDS 


CA-P GA-P 


TFN-P 


60 


220 


10 51 


0.9 


160 


10012 


2073 3131 
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tions. This constitutive model reproduces all the key fea- 
tures characteristic of these solutions with reasonable ac- 
curacy, including coil-stretch hysteresis (not shown here), 
and thus represents a significant improvement over con- 
stitutive equations such as the Oldroyd-B or FENE-P 
dumbbell equations typically used in the simulation of 
complex flows of dilute polymer solutions. 



VII. SUMMARY AND CONCLUSIONS 

The Gaussian Approximation was originally suggested 
for Hookean dumbbells and Rouse chains with HI by 
Ottingcr— and Wedgewood^, and has been shown to im- 
prove predictions by approximately accounting for the 
influence of fluctuations in HI. In this work, we use 
the same idea in the context of a free-draining bead- 
FENE spring model to first obtain a new approxima- 
tion for tackling the closure problem posed by the FENE 
nonlinearity. We show that the resulting equations in 
the approximate model can be interpreted in terms of a 
mean-field spring potential, and that this "FENE-PG" 
spring potential improves upon the well known FENE-P 
approximation in situations where spring force fluctua- 
tions are important. Furthermore, the use of additional 
diagonalization-of-normal-modes assumptions for simpli- 
fying the model equations lead to considerable savings in 
computation time, without significant loss in accuracy. 

The direct application of the Gaussian Approximation 
for bead-spring models simultaneously incorporating the 
effects of both FE and HI leads to an undesirable depen- 
dence of the predictions of equilibrium static properties 
on the hydrodynamic interaction parameter, h* . Instead, 
we use a two-step procedure to marry the FENE-PG 
mean-field spring potential with the Gaussian Approx- 
imation for HI to arrive at a closed set of equations for 
the second-moments, which is consistent with the equilib- 
rium Boltzmann probability distribution. This approxi- 
mate model accounts for the influence of fluctuations in 
both HI and the spring forces. 
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The method of combining the Gaussian approximation 
for HI with the FENE-PG approximation for HI sug- 
gested in this study is found to lead to predictions that 
compare well with the results of exact BD simulations 
for properties in steady and unsteady shear flows. The 
better agreement of the predictions of this approximate 
model — dubbed here as the "GA-PG" model — than those 
obtained with the FENE-P and/or Consistent- Averaging 
approximations emphasizes the importance of fluctua- 
tions in HI, and also highlights the role spring force fluc- 
tuations play at moderate and high shear rates. In ex- 
tensional flows, the GA-P model combining the Gaussian 
approximation for HI with FENE-P springs is found to 
perform better when chains are highly stretched, perhaps 
indicating that under such conditions, fluctuations in HI 
may still be important, even though local fluctuations in 
spring lengths are small. 

The predictions obtained with the Gaussian approxi- 
mation for HI in steady and unsteady shear and exten- 
sional flows confirms previous observations on the role 
of fluctuations in HI on macroscopic behaviour. When 
chains arc close to their equilibrium isotropic states, fluc- 
tuations in HI enhance the screening of the solvent veloc- 
ity gradient caused by HI. The role of fluctuations in HI 
reverses in situations where polymer chains experience 
significant stretching, and aid the penetration of the ve- 
locity field into the stretched polymer coil. Consequently, 
the solvent is able to engage a firmer "grip" on the chains 
and this tends to deform an already anisotropic chain fur- 
ther. Fluctuations in the FENE spring force arc found 
to increase the chains' resistance to stretching. 

When the stress predicted in BD simulations of FEBS 
chains during the start-up and following the cessation 
of a strong extensional fiow is plotted against predic- 
tions of the mean-squared end-to-end distance of the 
chains, or against the intrinsic birefringence, large hys- 
teresis loops were observed. The hysteresis loops pre- 
dicted with closure approximations were much smaller 
than those predicted with the simulations. Surprisingly, 
HI was found to have a negligible infiuence on the phe- 
nomenon of stress-conformational hysteresis observed in 
extensional flows. 

Closure approximations for FEBS chains with HI pre- 
dict coil-stretch hysteresis. In particular, predictions of 
hysteresis in the extensional viscosity by the GA-P model 
are in excellent agreement with the exact results obtained 
with BD simulations. It is observed that fluctuations in 
HI enhance the hysteretic effect, increasing the width of 
the coil-stretch hysteresis window on a plot of the steady- 
state extensional viscosity versus the strain-rate. 



by the expression. 



Hi-i-i. j — H, 



(Al) 



while the tensors Ay and AJ^ , which represent contri- 
butions due to fluctuations in HI, are given by. 



= ^'-J- ^^*^-r ■ ^"-^^ 



and, 



r,s — 1 



A^= ^(i..a:j:qt, 



(A2) 



with, 

(A3) 

The second-rank tensorial function ffj/^ = JT(S'*^) is 
related to the average of the HI tensor through. 



(Cn,^) = V2h*H{s*) , 



(A4) 



whereas the fourth-rank tensorial function Kf^i^ 
K{S*^^) is related to the HI tensors through. 



3^/2/1* 



In the equations above. 



(A5) 



(A6) 



Here, 



1, if yU < z < , 

Ae;,, = <|-i, ifz.<z<^, (A7) 

, otherwise, 

is the onc-dimensional "box" function. This function can 
be succinctly expressed in terms of the Heaviside step 
function. 



^ / N 1, if n > m / . „s 

9(771, n) = <^ ' - (A8) 

otherwise , 



APPENDIX A: Tensors associated with hydrodynamic 
interactions 



as Aej,^ = 6(^,i) - e(7/,i) = -AOj,^. Further, the 
dyad 



The tensor Ajj in Eq. H48|l . which describes contri- 
butions arising from averaging the Osccn-Burgers tensor 
with the non-equilibrium distribution function, is given 



max(^',^t) — 1 A^^ 

J2 Q^Q^ = E Ae;f,Q,Q, , (A9) 

i,j — min(L',/i,) i,j —1 
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where the two-dimensional box function, 



(AlO) 



This function is symmetric with respect to an exchange 
of the indices fi and v, and/or i and j. 

When the Oseen-Burgers form of the HI tensor is used, 
it can be shown that^iSS 



H, 



2(2' 



(All) 



and 
K 



(27r)3/2 J 



rfc 5 



fcfc\ 1 

— j fc exp(- - fefe : S*^) dk , 
(A12) 



Thus, the functions H^^ and Kufn and therefore Aij and 
Ay are completely expressible in terms of the second 
moments of the Gaussian distribution. 

The inner product of K with a second-rank tensor a 
is defined in component form as 



7,(5=1 



(A13) 



and, 



(a: K)ai3 = ^ as^K^Sa/i for a, (3 = 1,2,3 . (A14) 

7,(5=1 

Further, since the inner product 
1 / kk\ , 



for all k, we have 



K: 8 = 0. 



(A15) 



A direct consequence of this result is that the A ten- 
sors defined through Eqs. I|A2|I and IjASp . vanish under 
isotropic equilibrium conditions. 
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